
%% 测试参数设置
m = 1000;
n = 120;
r = 100; 
A = randn(m, r) * randn(r, n); 
T = eye(m); 
%% 算法执行与计时
fprintf('===== 算法性能对比 =====\n');

% MATLAB内置pinv
tic;
A_pinv = pinv(A);
time_pinv = toc;

% 自定义QRIL算法
tic;
Beta = QRIL_inverse.QR_schmidt(A, T);
time_qril = toc;

% 计算误差
error_pinv = norm(A*A_pinv - T, 'fro');
error_qril = norm(A*Beta - T, 'fro');
diff_pinv = norm(Beta - A_pinv, 'fro')/norm(A_pinv, 'fro');

%% 结果输出
fprintf('矩阵大小: %d×%d (秩=%d)\n', m, n, rank(A));
fprintf('---------------------------------\n');
fprintf('方法\t\t时间(ms)\t误差(Fro)\t与pinv差异\n');
fprintf('---------------------------------\n');
fprintf('pinv\t\t%.2f\t\t%.2e\t基准\n', time_pinv*1000, error_pinv);
fprintf('QRIL\t\t%.2f\t\t%.2e\t%.2e\n', time_qril*1000, error_qril, diff_pinv);

%% 数值稳定性测试（条件数影响）
fprintf('\n===== 数值稳定性测试 =====\n');
cond_numbers = [1e2, 1e4, 1e6, 1e8, 1e10, 1e12];
results = zeros(length(cond_numbers), 3);

for i = 1:length(cond_numbers)
    cond_num = cond_numbers(i);
    
    % 正确构造指定条件数的矩阵
    [U, ~, V] = svd(randn(m, n), 'econ');
    s = linspace(1, 1/cond_num, min(m,n));
    A_cond = U * diag(s) * V';
    
    % 计算伪逆
    A_pinv_cond = pinv(A_cond);
    Beta_cond = QRIL_inverse.QR_schmidt(A_cond, eye(m));
    
    % 记录误差
    error_pinv_cond = norm(A_cond*A_pinv_cond - eye(m), 'fro');
    error_qril_cond = norm(A_cond*Beta_cond - eye(m), 'fro');
    diff_cond = norm(A_pinv_cond - Beta_cond, 'fro')/norm(A_pinv_cond, 'fro');
    
    results(i, :) = [error_pinv_cond, error_qril_cond, diff_cond];
    fprintf('条件数=%.0e: pinv误差=%.2e, QRIL误差=%.2e, 相对差异=%.2e\n', ...
            cond_num, error_pinv_cond, error_qril_cond, diff_cond);
end

%% 可视化结果
figure;

% 误差比较
subplot(1, 2, 1);
semilogy(cond_numbers, results(:, 1), 'b-o', 'LineWidth', 1.5, 'DisplayName', 'pinv误差');
hold on;
semilogy(cond_numbers, results(:, 2), 'r-s', 'LineWidth', 1.5, 'DisplayName', 'QRIL误差');
xlabel('条件数');
ylabel('误差 (Frobenius范数)');
title('伪逆计算误差比较');
legend('Location', 'northwest');
grid on;
set(gca, 'XScale', 'log');
set(gca, 'YScale', 'log');

% 算法差异
subplot(1, 2, 2);
semilogy(cond_numbers, results(:, 3), 'm-d', 'LineWidth', 1.5);
xlabel('条件数');
ylabel('相对差异 (Frobenius范数)');
title('与pinv结果的相对差异');
grid on;
set(gca, 'XScale', 'log');
set(gca, 'YScale', 'log');

%% 验证Moore-Penrose条件
fprintf('\n===== Moore-Penrose条件验证 =====\n');
test_matrices = {A, randn(m, n)}; % 原始矩阵和随机矩阵

for i = 1:length(test_matrices)
    A_test = test_matrices{i};
    
    % 计算伪逆
    A_pinv_test = pinv(A_test);
    Beta_test = QRIL_inverse.QR_schmidt(A_test, eye(m));
    
    % 验证四个Moore-Penrose条件
    mp1_pinv = norm(A_test*A_pinv_test*A_test - A_test, 'fro');
    mp1_qril = norm(A_test*Beta_test*A_test - A_test, 'fro');
    
    mp2_pinv = norm(A_pinv_test*A_test*A_pinv_test - A_pinv_test, 'fro');
    mp2_qril = norm(Beta_test*A_test*Beta_test - Beta_test, 'fro');
    
    mp3_pinv = norm((A_test*A_pinv_test)' - A_test*A_pinv_test, 'fro');
    mp3_qril = norm((A_test*Beta_test)' - A_test*Beta_test, 'fro');
    
    mp4_pinv = norm((A_pinv_test*A_test)' - A_pinv_test*A_test, 'fro');
    mp4_qril = norm((Beta_test*A_test)' - Beta_test*A_test, 'fro');
    
    fprintf('\n矩阵 %d:\n', i);
    fprintf('条件1: pinv=%.2e, QRIL=%.2e\n', mp1_pinv, mp1_qril);
    fprintf('条件2: pinv=%.2e, QRIL=%.2e\n', mp2_pinv, mp2_qril);
    fprintf('条件3: pinv=%.2e, QRIL=%.2e\n', mp3_pinv, mp3_qril);
    fprintf('条件4: pinv=%.2e, QRIL=%.2e\n', mp4_pinv, mp4_qril);
end


% m = 100;
% n = 100;
% k = 30;
% A = rand(m, n);
% beta_std = rand(n, k);
% T = A * beta_std;
% 
% % 计算 QR_schmidt 的时间和结果
% tic;
% Beta = QRIL_inverse.QR_schmidt(A, T);
% elapsed_time = toc;
% time_QR_sch = elapsed_time;
% norm_QR_sch = norm(A * Beta - T, "fro")


% m = 100;
% n = 120;
% A = rand(m, n);
% d = min (m,n)
% T = diag(d);
% Beta = QRIL_inverse.QR_schmidt(A, T);
% norm(A * Beta - T, "fro")
% beta_std = rand(n, k);
% 
% T_true = A * beta_std;


% if m>n
%     Beta = QRIL_inverse.QR_schmidt(A', T);
%     Beta = Beta';
% else
%     Beta = QRIL_inverse.QR_schmidt(A, T);
% end



% tic;
% if m>n
% 
% Beta = QR_schmidt(A, T);
% 
% else
% 
% end
% 
% elapsed_time = toc;